Volume 108, Number 6, November-December 2003 

Journal of Research of the National Institute of Standards and Technology 

[J. Res. Natl. Inst. Stand. Technol. 108, 413-427 (2003)] 

Simulation of an Austenite-Twinned-Martensite 

Interface 



Volume 108 



Number 6 



November-December 2003 



A.J. Kearsley and L. A. Melara 
Jr. 

National Institute of Standards 
and Technology, 
Gaithersburg, MD 20899-0001 

ajk@nist.gov 
luism@nist.gov 



Developing numerical methods for pre- 
dicting microstructure in materials is a 
large and important research area. Two 
examples of material microstructures are 
Austenite and Martensite. Austenite is a 
microscopic phase with simple crystallo- 
graphic structure while Martensite is one 
with a more complex structure. One 
important task in materials science is the 
development of numerical procedures 
which accurately predict microstructures 
in Martensite. In this paper we present a 
method for simulating material microstruc- 
ture close to an Austenite-Martensite inter- 
face. The method combines a quasi- 
Newton optimization algorithm and a non- 
conforming finite element scheme that 



successfully minimizes an approximation 
to the total stored energy near the interface 
of interest. Preliminary results suggest that 
the minimizers of this energy functional 
located by the developed numerical algo- 
rithm appear to display the desired charac- 
teristics. 
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1. Introduction 

Simulation of microstructures is important for deter- 
mining the behavior of complex materials [1,2,3,4,5,6]. 
Applying macroscopic loads to some materials can 
reveal their mechanical properties. These properties are 
associated with the microscopic structures created due 
to the material deformation. Microscopic structures of 
interest are often found near internal surfaces which, in 
turn are found in systems of liquid-liquid, solid-liquid 
and solid-solid boundaries [7]. 

An Austenite-twinned-Martensite interface is an 
example of an internal surface. This interface conjoins 
two phases: Austenite and Martensite. Figure 1 shows 
an example of the twinned-Martensitic microstructures. 
These microscopic structures are the result of thermal 
or mechanical loading of the material. When a material 
is deformed, the structural organization of the atoms is 
rearranged. At the microscopic level, this rearrange- 



ment can create the crystallographic patterns shown in 
Figure 1, for example. These patterns represent 3- 
dimensional microscopic structures. There are many 
different ways to rearrange the atoms that make up a 
material and each variant corresponds to a particular 
arrangement. We focus on the two specific variants of 
Martensite distinguished by the alternating shaded 
areas in Figure 1 . 

One important mechanical property associated with 
an Austenite-twinned-Martensite interface is the Shape 
Memory Effect (SME) characteristic of Shape Memory 
Alloys (SMAs). SMAs are materials that, when 
deformed for an indefinite period of time, return to their 
original shape. The deformed state is a metastable state 
[8]. In metastability, the total stored energy is at a local 
minimum, thus requiring energy to induce a transfor- 
mation. In SMAs, heat will trigger the SME which in 
turn, returns the SMA to its original shape. An 
Austenite-twinned-Martensite interface is observed in 
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Fig. 1. Twinned-Martensite is shown to the lower left of the diagonal line. The field of view 
is in the order of |am. 



materials when they are in a metastable state. 
Therefore, to simulate Martensite, we minimize the 
total stored energy near an internal surface. In this 
paper, we present a numerical technique used to obtain 
a solution which represents the spatial structure of the 
twinned Martensite. The technique employs a Q^i finite 
element discretization of a function representing total 
energy and a limited memory quasi-Newton method to 
minimize the resulting approximate total energy func- 
tion. The gradient of the function is calculated by 
computing the partial derivatives using finite differ- 
ences. The result is an apparently robust method for 
approximating what are usually difficult minimizers to 
locate. 

1.1 Total Stored Energy 

A myriad of research methods can be found which 
focus on the simulation of twinned Martensite, and 
related microstructures, e.g. Refs. [9,10,11,12,13,14]). 
Previous numerical work includes minimization of the 
total stored energy using the conjugate gradient meth- 
ods [15], the method of steepest descent [15], and a 
descent method [16]. Further references using and 
studying the use of finite element methods to simulate 
Martensitic micro structures can be found in 
[17,18,19,20,21,22,23,24]. We are interested in simu- 
lating the microstructures in the metastable state, thus 
our work will focus on the static problem. 



The functional representing the total stored energy is 
taken from work by Kohn-Muller in [10,11]. The 
Martensite region is represented by the domain 
n=(0,L)x(0,K), Letx = (x,j)e I3czM',w:M'^M 
and u = u{x,y). The double-well energy function is 
given by 

j(u)=i(d^ur+((d^ur-]j<ix+€jjd^u\<ix (1) 

where u equals at x = 0. The x = boundary corre- 
sponds to the internal surface. The function w = at 
X = 0, due to the constraint of elastic compatibility 
[ 1 0, 1 1 ] . The first integral is the elastic energy and the 
last integral is the surface energy. We seek a function 
u e yV^'"^ (13) which minimizes Eq. (1) where 

and v = atx = 0.} 

1.2 Elastic Energy 

The 2-D scalar model of elastic energy, in Eq. (1) is 
l[(d^uf+{(d,Mf-lj]dx (2) 

This term is minimized by a function u such that 



VW: 



±1 



(3) 
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Each of the gradients in Eq. (3) corresponds to a stress- 
free state in one of the two distinct variants of 
Martensite [10,11]. Examples of functions satisfying 
Eq. (3) are plotted in Fig. 2. 

In Fig. 2 functions are piecewise linear in thej^-direc- 
tion and constant along the x-direction. Requiring that 
u equal zero at x = creates more oscillations. This is 
similar to the behavior of Martensite near an internal 
surface. The function with Amplitude 1 is closest to 
zero, on average, than those with Amplitudes 2 and 3. 

The number of oscillations is directly related to the 
number of discontinuities in dyU. We note that the func- 
tion with Amplitude 1 also has the largest number of 
discontinuities, occurring at the peaks of u, where dy u 
changes from +1 to -1 (or vice-versa). Therefore, a 
desirable minimizer u has the characteristics of the 
function with Amplitude 1 . The consequence of u satis- 
fying Eq. (2) and w = at x = is the occurrence of 
more discontinuities in dyU. At the minimizer u, howev- 
er, the surface energy penalizes these discontinuities. 

1.3 Surface Energy 

We employ the surface energy presented by Kohn 
and Miiller in Ref [10,11], 



£] \du\6x. 



(4) 



where 8 is a constant. 

The definition proposed by Kohn and Miiller 
replaces the integral of Eq. (4) in the j^-direction with 
the total variation of dyU over (0,^ [10,11]. This 
change designates the role of the surface energy at 
the minimizer w as a "counter" because ^Jf I 3^w | dy 
counts the number of discontinuities in dyU, where 
dyU^ {±1}- Thus, at the minimizer u of Eq. (1), the 
surface energy Eq. (4) exhibits opposite behavior to Eq. 
(2) since we are minimizing Eq. (1). To illustrate this 
counting role of Eq. (4), we briefly review functions of 
bounded variation. 

1.4 Functions of Bounded Variation 

The definition of functions of bounded variation is 
taken from Refs. [25,26,27]. Let 

be any partition of (0, K) and let 

|y|=max|3;.^,-3;.|. 
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Fig. 2. Oscillatory behavior of minimizer u. 
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The partitioning T is a collection of points yp 7 = 0, 1 , 
. . . , /, such that yQ = and j/ = K. For each partitioning 
y, let 

Sy(x) = ^1 dyU(x, yj^,)-dyU(x, J.) I, (5) 

7=0 

for a fixed xe (0, Z). The total variation of dyU over 
(0, K) is 



supSy(x), 



(6) 



where the supremum is taken over all countable 
partitions T of (O^K). If < Sy(x)<+oo, then 0< 
^1^2 1 3^^ I dx<+oo. If ej^l dyyU I dx<+oo, then3^w is 
a function of bounded variation. 

Next, we describe the substitution of Jf | 9^w | dj 
by Eq. (6). We assume that 3^w e C\Q). Using Eq. (5) 
and the Mean Value Theorem, we have for some 
Cye (ypyjH%J= 1, 2, ..., /that 
/-I 



7=0 

/-I 



^Xl^>T^(-^'^y)l(->'7+i"-^7)' 



(7) 
(8) 



7=0 



and we obtain by Theorem 2.9 in Ref [25]: 



I— I 
supSy (x) = lim Sy (x) = lim Xl a^w(x, g^) | (^^^^^ -3;^) 

= JolMdj (9) 

The assumption that u e C\Q.) is stronger than the 
assumtion that u e ^V^'"^ (Q) ; however, the equivalence 
in Eq. (9) is possible since the changes in dyU do not 
occur in subintervals of measure zero, [9,10,11,28]. 
Figure 3 shows that ^J^ | d^u \ dy counts the number 
of sign changes in dyU. Let (0, K) = (0, 1). 

In Fig. 3a, the piecewise linear function uQc,y) is 
plotted over (0, 1). In Fig. 3b, the function dyuQc,y) is 
plotted over (0, 1). Clearly, the number of discontinu- 
ities in 3^w is 5 from Fig. 3b. For a discontinuity at the 
partition point >?y, we adopt the following convention: 

w(x, v)-w(x, V.) 
d,uix,yj)^ lim -^^ (10) 

y-^y,.y>y, y-yj 

Now consider (y,_i, j,+i) = O^i, Jy)Ufe J'y^i) i" Fig' 3b. 
Evaluating at the endpoint of the two subintervals gives 

\d^u{x,yj)-d^u{x,yj_{)\=\\-{-\)\ =2 and (11) 
\^Ax,yj^,)-d^u{x,yj)\= 1 1-1 1 =0. (12) 




(b) 







-1 -- 



Fig, 3, Surface energy term as a counter. Pictorial example illustrat- 
ing the role of surface energy \\ \ dyyU | dy as a counter term. The 
notation u{-, y) means that the graph shows the >'-dependence only. 



A jump occurred in (j^^i, jy) because dyU is discontinu- 
ous at a point in this subinterval. Based on our conven- 
tional definition of 3^ w(x, yj) at a point 3^^ in Eq. (10), we 
say a "jump" occurs at;^^ if 3^w is discontinuous at this 
point. This "jump" refers to the change in values of 3^w 
aXyj from +1 to -1 or vice-versa. 

This jump gives the nonzero value in Eq. (1 1). In the 
subinterval (ypyj^i), no jump occurs in dyU, therefore 
the difference in Eq. (12) is zero. We compute Eq. (9): 

\\\dyyU{x,y)\dy = 2x5 =\Q, 

where the magnitude of the jump is 2 and the number 
of discontinuities is 5. Table 1 shows the number of dis- 
continuities for the functions presented in Fig. 2. 

Table 1, Number of "jumps" in d u 



Function 



Jo I dyyU I d>' 



Discontinuities 



Amplitude 1 
Amplitude 2 
Amplitude 3 



24 

12 

6 



12 
6 
3 
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Therefore, we have: 

J I dy^,u(x,y) \dy=2x (number of discontinuities of 3^,ii), 

and the surface energy succeeds in serving as a "count- 
er." The next section describes the numerical approxi- 
mation to Eq. (1). 



2, Description of Method 

The total stored energy Eq. (1) consists of elastic 
plus surface energy. First, we describe the implementa- 
tion of the elastic energy followed by a discussion of 
the surface energy implementation. We computed Eq. 
(2) using affine finite elements, [29]. We chose the 
finite element space: 



£i =span{l,x,>',XF}, 



(13) 




because functions in this polynomial space lie in 
^^'"^(Q), which satisfy the criteria in Eq. (3) for a min- 
imizer u of Eq. (1). The elements are rectangles with 
the function values at the vertices as degrees of free- 
dom. The parent element is Q = (0, 1) x (0, 1) and for 
some A: e N, the reference element is 

P 

where h^ is the mesh size along the x-axis and h2 is the 
mesh size along thej^-axis. The degrees of freedom of 
the elements are the function values at the vertices, and 
for Q we label them w(ai), w(a2), u{^^, and u{^^. Figure 
4 shows the parent element with the vertices also denot- 
ed by concentric circles. 

2.1 Affine Finite Elements 



Fig. 4. Parent Element Q with four degrees of freedom labeled w(a|), 
w(a2), w(a3), and u{2i^. 



where x = x(x, y), y = y(x, y). The nth basis function in 
Qis 

and it satisfies the property that 0„(aj = 6^^, for n,m = 
1, 2, 3, 4, where 5^^ is the Kronecker delta function. 

To compute the gradient of u and to integrate u, we 
need: VF, detVFand VF~^. The gradient of Fis 



VF = 



h, 
h. 



(16) 



(17) 



Let F : Q^Qk be given by 



>' 





V 



0\x 



y 






(14) 



We simplify, by letting F = F{x,y). We use the affme 
map Eq. (14) to evaluate the function w, to compute Vw 
and to integrate. The function u is spanned by the four 
basis functions defined in each element g^. Using Eq. 
(14), we can determine all basis functions by the four 
basis functions in Q. Thus, the approximation to u is: 

w(x,3;) = Ui<^i(x,j)+U202(-i,j)+U3<^(^j)+U4 04(^3;), 

(15) 



and, we have 

Using Eqs. (15), (16), and (17) we have 

Vm(x,j) = VF-'£u„V0„(x,>'). 



To integrate the elastic energy over Qi, we use Eqs. 
(15)-(17) and obtain 



J^^ {d,uf6xAy = J j ^-jj^At Cx, y) det VFd5d>, 

=ji|l»-.5A(^'i^) j hih^^y- (18) 
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The procedure is the same for 

2.2 Computation of Surface Energy 

We evaluate the surface energy term in Eq. (4) using 
Eq. (6), approximating d^u using backward finite differ- 
ences, 



dJ(u) , 



and 



^yu(x,yj)' 



^yti(x,yj^,) 



^(x,yj)-^(x.yj-i) 



^(x.yj^i)-^(x.yj) 



Here, 



JJ9^w|d>'= X, \dyU(x,yj^,)-dyU(x,yj)\ (19) 



where N(h2) is the number of nodes along the j^-direc- 
tion. We can integrate along the x-direction to complete 
the approximation of Eq. (4). 



Giu)=-Jiu+t<pn__,=i-yv<pdx, 

nf •'^^ Hv;y 



d^ 



dVu 



-\ div (pdx. 



where the gradient is div (dJ/dVu). Let J(u) be the dis- 
cretized representation of the total energy Eq. (1) for 
u e M^^^^, where N(h) is the total number of nodes in Q. 

With this numerical approximation to VJ(u), we use 
the limited memory BFGS method [32] to numerically 
minimize J(u). This inexpensive quasi-Newton method 
seeks to build an approximation H^^^ to the inverse of 
the Hessian matrix of second derivatives of J at the 
point u^^l These inexpensive approximations are con- 
structed from a small number of vectors updated at 
every iteration. 

The k+ I step of the limited memory BFGS method 
is given by: 

where the parameter d^^ is a line search parameter and 
is computed at each iteration using a line-search proce- 
dure. 

The Hessian matrix B^^^ satisfies the secant equation: 



B(%('^=/'\ 



where 



3. Numerical Implementation 

The computation of the integral in Eq. (18) is done 
exactly in the parent element Q. The symbolic package 
Maple^ [30] was employed to integrate this term exact- 
ly. The algebraic expression consisting of nodes Uj and 
mesh sizes hy and h2, is evaluated over each element. 

The gradient of Eq. (1) is approximated using cell- 
centered finite differences. Let w Gly^''^(Q). The 
Gateaux derivative is: 

G(u) = \im-[j(u + t(p)-J(u)] (20) 

where <p e "H^^'"^ (Q) [31]. The linear operation G is the 
directional derivative of Jin the direction <p. From the 
Euler-Lagrange equations, we see that the gradient of J 
is difficult to compute. Considering the Euler-Lagrange 
equation we have. 



^ Certain commercial equipment, instruments, or materials are iden- 
tified in this paper to foster understanding. Such identification does 
not imply recommendation or endorsement by the National Institute 
of Standards and Technology, nor does it imply that the materials or 
equipment identified are necessarily the best available for the pur- 
pose. 



Ak) _ (,Xk) 



= (w^^^-w^^-^0, and 
yW^VJ(w^'^)-VJ(w^'-^^). 

We employ a limited memory BFGS method because 
the Hessian matrix B^^^ is expensive to compute since 
we are solving a large scale optimization problem 
[32,33]. This is done by storing a certain number of 
vector pairs {s^^\ y^^}. At each new iterate, the oldest 
vector pair is deleted and replaced by the new vector 
pair, thus preserving curvature information only from 
the most recent iterations. The formula for updating the 
inverse of the Hessian matrix H^^^ = (B^^^)'^ is: 

Jj(k+l) _/y(k-m) yik-l) sT r fj{k) \0 ry{k- m) y{k-\) ^ 

j^ {k-m),Y{k-m+\) y{k-\)Y ^{k-m),^{k-m)-yT ,y{k-n^\) yik-l) \ 

^ {k-m+\),y{k-m + 2) y{k-\)-yT ^{k-m^\),^k-m+\)-yT ,y{k-ni^l) ^^1)\ 

+p(^-l)g(^-l)(g(^-l))^_ 
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where (If^^f is some initial Hessian approximation, 

and 

1 



(f^)\T Jk) • 



i/'^S 



4. Results 



One final and curious observation is the appearance 
of a "diagonal band" structure in Figs. 5-10. It would 
seem that this band, which is apparent for solutions on 
various grid sizes, is related to the competing roles of 
the surface and elastic energies is and may be associat- 
ed to the "equipartitioning of energy" principle pro- 
posed in Kohn and Miiller in Ref [10,11]- Future 
research plans include a deeper investigation of this 
diagonal-band structure. 



In this section, results computed on the various grids, 
30 X 30, 60 X 60, and 120 x 120, are presented. In all 
minimization, the initial iterate u^ is a small (~10~^) per- 
turbation of pure Austenite, which we take to be w = 0. 
Density plots of the minimizer u are shown in Figs. 5- 
10 and demonstrate the existence of the desired 
microstructures in the twinned-Martensite phase. These 
figures also demonstrate the role of the surface energy 
as a penalization term of Eq. (1). We seem to be able to 
control the number of microstructures by varying the 
values of £, with a larger value of £ resulting in a small- 
er number of microstructures in Martensite and vice- 
versa. Figures 11, 12, and 13 show plots of the profiles 
of minimizer w atx = 4. Figures 11, 12, and 13 exhibit 
a larger number of discontinuities in dyU for £=2x 
10^ than for £=2x 10"^. In Figs. 5-9 one sees the 
effect of grid-size comparing results on 60 x 60 grid 
with 120 X 120 grid for £=2x 10^ through £=2x 
10"^. The black and white stripes correspond directly to 
the saw-toothed behavior of w at x = 4. For small £ val- 
ues, the number of discontinuities is close to the num- 
ber of partitions along the j^-direction. Tables 2-4 show 
the energy values for the minimizers shown in Figs. 5- 
10. The column with E(u) corresponds to the elastic 
energy values of the minimizer while the column with 
S(u) corresponds to the surface energy values of the 
minimizer. The energy values presented in Tables 3 and 
4 appear consistent with these conclusions. The larger 
the £, the larger the value of the surface energy at the 
minimizer. In Table 3 the value of the surface energy 
increases by a larger order than that of the elastic ener- 
gy 



Table 2, Energy values for various £ on 120 x 120 grid 

£ E{u) £S{u) Total energy 



2x10^ 
2 X 10"^ 
2x10^ 
2 X 10"^ 
2x10"^ 



6.3014 


0.569 le-2 


6.1122 


0.5698e-l 


7.0859 


0.5474 


8.2222 


4.5200 


23.3218 


1.2559 



6.3071 

6.1692 

7.6334 

12.7423 

24.5770 



Table 3, Energy values for various £ on 60 x 60 grid 

£ E{u) £S{u) Total energy 



2x10^ 


4.3684 


0.3120e-2 


4.3716 


2 X 10"^ 


6.3005 


0.2808e-l 


6.3286 


2x10^ 


6.3010 


0.2798 


6.5809 


2 X 10"^ 


5.7738 


2.7356 


8.5095 


2x10"^ 


3.7945 


14.6825 


18.477 



Table 4, Energy values for various £ on 30 x 30 grid 



£ 


E(u) 


£S(U) 


Total energy 


2x10^ 


4.6093 


0.0015 


4.6108 


2 X 10"^ 


4.9832 


0.0145 


4.9977 


2x10^ 


6.1073 


0.1326 


6.1499 


2 X 10"^ 


5.5335 


1.3246 


6.8581 


2x10"^ 


8.0650 


9.6770 


17.7420 
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Fig. 5. Comparison of density plots of w for values of £ = 2 x 10 on 60 x 60 grid (bottom) vs 120 
X 120 grid (top). 



£ = 2.00d-5 




e=2.00d-5 




Fig. 6. Comparison of density plots of w for values of £ = 2 x 10 on 60 x 60 grid (top) vs 120 x 120 
grid (bottom). 



420 



Volume 108, Number 6, November-December 2003 

Journal of Research of the National Institute of Standards and Technology 



e = 2.00d-4 
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16 



£ = 2.00ci-4 




Fig. 7. Comparison of density plots of w for values of £ = 2 x 10 on 60 x 60 grid (top) vs 120 x 120 
grid (bottom). 



e=2.00cl-3 




£ = 2.00ci-3 




Fig. 8. Comparison of density plots of u for values of £ = 2 x 1 ^ on 60 x 60 grid (top) vs 1 20 x 1 20 
grid (bottom). 
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£ = 2.00d-2 




£ = 2.00d-2 




Fig. 9. Comparison of density plots of u for values of £ = 2 x 1 on 60 x 60 grid (top) vs 1 20 x 120 
grid (bottom). 
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Fig. 10. Density plot of solution w for £ = 2 x 10"^, ..., 2 x 10 ^ on 30 x 30 grid. 
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Fig. 11. Comparison of w(4,>') for various values of £= 2 x 10 , ..., 2 x 10 on 120 x 120 (left column) grid vs 60 x 60 grid (right column). 
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t:=2.00cJ-3 



E - 2.00d-3 





t:-2.00d-2 



i:\ - 2.00d-2 



0.5 1 



25 3 35 
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Fig. 12. Comparison of u{4, y) for various e = 2 x 10 and e= 2 x 10 on 120 x 120 grid (left column) vs 60 x 60 grid (right column). 
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Fig. 13. Profiles of w(4,>^) for £ = 2 x 10 , ..., 2 x 10 on a 30 x 30 grid. 
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